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Understanding relationship between noisy dynamics and biological network architecture is 
a fundamentally important question, particularly in order to elucidate how cells encode and 
process information. We analytically and numerically investigate general network architec- 
tural conditions that are necessary to generate stochastic amplified and coherent oscillations. 
We enumerate all possible topologies of coupled negative feedbacks in the underlying bio- 
chemical networks with three components, negative feedback loops, and mass action kinetics. 
Using the linear noise approximation to analytically obtain the time-dependent solution of 
the master equation and derive the algebraic expression of power spectra, we find that (a) 
all networks with coupled negative feedbacks are capable of generating stochastic amplified 
and coherent oscillations; (b) networks with a single negative feedback are better stochas- 
tic amplified and coherent oscillators than those with multiple coupled negative feedbacks; 
(c) multiple timescale difference among the kinetic rate constants is required for stochastic 
amplified and coherent oscillations. 



I. INTRODUCTION 



Many biological processes such as gene regulation, cellular signaling events, and metabolic 
reactions typically involve a large number of heterogeneous biochemical constituents and their 
random interactions at and across multiple temporal and physical scales [U [2] . They can be mapped 
onto complex networks where a pair of constituents are linked if they interact with each other 
physically or chemically [3^ 4J. The topological properties of the large-scale biological networks 
have been extensively investigated, yet, the dynamic behaviors of and on such large-scale networks 
remains unexplored and of great challenge despite its importance in understanding many biological 
timing mechanisms such as cell cycle, circadian rhythms, neural rhythms, cardiac rhythms, and 
hormone rhythms [5, 6J. The technical challenges are primarily due to not only a lack of information 
of accurate network structure and kinetics, but also a sheer number of biochemical components 
and their random interactions that requires a new mathematics that can handle high-dimensional 
stochastic nonlinear dynamical systems. Our paper is to provide an insight into such an important 
stochastic nonlinear dynamical behaviors of small-size biological networks, namely coupled network 
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motifs. 

Network motifs are recurring subgraphs whose appearance in certain biological networks is much 
more frequent than would be expected in random networks and can be considered as simple building 
blocks from which the network is built 0^9]. There has been extensive research on identification 
of network motifs and their dynamical and functional roles in biological networks [8Hl3j. Most of 
the previous studies have been limited to investigation of the dynamics of basic network motifs in 
isolation. In real biological networks, those individual network motifs are rather coupled with each 
other and embedded in a larger network. Therefore, it is crucial and even imperative to study the 
dynamical behaviors of coupled network motifs. 

In Ref . [HI [15] , Tyson et. al. classified the biochemical oscillators according to the topology 
of coupled negative and positive feedback loops in the underlying regulatory networks with three- 
components and discussed the general network architectural requirements for oscillations: negative 
feedback and time-delay by a series of intermediate components. In this regard, a network with 
two-component negative feedback loop without explicit time-delay cannot generate oscillation, 
but a two-component negative feedback loop coupled with a positive feedback loop is able to 
generate oscillation. Here positive feedback plays a role of time-delay. Tsai et. al. uncovered an 
additional role of a three-component network with interlinked positive and negative feedback loops: 
frequency tunable oscillator pO]. Also, it was shown that the dynamical role of coupled negative 
feedbacks is to realize enhanced homeostasis whereas coupled positive and negative feedbacks 
properly modulate signal responses and effectively deal with noise |16fll9j . The previous dynamical 
studies on coupled feedbacks neglected stochastic fluctuations that are prevalent and sometimes 
dominant in intracellular biochemical processes [ 211 122j . 

Stochastic fluctuations can enhance the net order and induce stochastic resonance [23] where 
white noise amplifies a weak periodic external signal and the signal to noise ratio is maximized at 
the optimal noise strength. The dynamical systems that are not subject to any periodic external 
forcing can exhibit the similar phenomena coined as "coherence resonance" or "noise induced os- 
cillations" |24H27| [3TVl33j . In the latter case, intrinsic noise can induce stochastic amplified and 
coherent oscillations out of a stable dynamical system. Many excitable systems are known to exhibit 
noise-induced oscillations |26| [28H30j and this excitability is mistakenly assumed to be necessary 
for noise-induced oscillations. However, even in non-excitable systems such as a predator-prey 
model, particularly where the deterministic dynamical systems are stable in the entire parameter 
space, noise give rise to stochastic oscillations with high amplification and coherence \27\ \31 \ [32]. 
Presently, there exist no theoretical work about the rigorous conditions for noise-induced oscilla- 
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tions, especially when the deterministic dynamical systems have no center and thus no limit cycle 
in the entire parameter space. 

In this paper, we analytically and numerically investigate general network architectural con- 
ditions that are necessary to generate stochastic amplified and coherent oscillations. We classify 
stochastic oscillators according to the topology of coupled negative feedback loops in the underlying 
biochemical networks with three components and furthermore identify the best network designs 
that can generate stochastic amplified and coherent oscillation. In the following sections, we intro- 
duce an exhaustive list of nine networks with three components, negative feedback loops, and mass 
action kinetics. We use the linear noise approximation to analytically obtain the time-dependent 
solution of the master equations and derive the algebraic expression of power spectra. The signal 
to noise ratio measured from a power spectrum is used to determine if the particular topology of 
interconnected negative feedbacks is a good design for stochastic oscillation. We find that (1) all 
nine networks with three-component negative feedbacks and with mass action kinetics are capable 
of generating stochastic amplified and coherent oscillations; (2) networks with a single negative 
feedback are better stochastic amplified and coherent oscillators than those with multiple coupled 
negative feedbacks; (3) multiple timescale difference among the kinetic rate constants is required 
for stochastic amplified and coherent oscillations. 

II. METHODS 
A. Models 

We consider simple and idealized cell signaling networks consisting of 3 nodes {X, Y, Z) . The 
nodes are the chemical species participating in the networks. The chemical reaction occurring 
among the nodes is represented by an edge connecting them. The reaction can be either an 
activation (e.g., X ^Y) or an inhibition (e.g., X H Y). 

In this work, we pay a particular attention to the networks that are entirely composed of the 
negative feedback loops (NFBLs) . The simplest networks among all are the ones with 3 nodes and 
3 edges. There are only two possible configurations, i.e., i) NFBL network 1: X^Y^Z-\Xox 
ii) NFBL network 2: X-\Y-\Z-\Xso that they completes two different types of 3 component 
single NFBLs. If we consider more edges, we can construct multiple feedback loop networks. All 
the multiple NFBL networks can be built either based on NFBL network 1 or 2. We show the 
names and network configurations of all the networks in Table |I] and its graphical representation 
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in Fig. [T] They constitute the exhaustive hst of the distinctive networks consisting of the negative 
feedback loops alone with 3 nodes. Among them, 6 multiple NFBL networks (NFBL networks 3, 4, 
5, 7, 8, and 9) are based on NFBL network 1. In particular, NFBL networks 3, 4, and 5 consist of 
double NFBLs and NFBL networks 7, 8, and 9 triple NFBLs. There are 2 multiple loop networks 
(NFBL networks 6 and 10) based on NFBL network 2. NFBL networks 6 and 10 are made of 
double and triple NFBLs respectively. 




FIG. 1: Exhaustive list of all distinct networks with 3 components and negative feedback loops. 

The networks are numerically labeled and ordered by the number of negative feedback loops: 
single NFBL in the top row, two NFBLs in the middle row, and three NFBLs in the bottom row. 
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TABLE I: Network Names and their topology 



Names of network 


network topology 


NFBLl 


X ^ Z ^ X 


NFBL2 


X ^ Z ^ X 


NFBL3 


X^Y^Z^XanAY^X 


NFBL4 


X^Y^Z-\X and Z-\Y 


NFBL5 


X^Y^Z^X and X^Z 


NFBL6 


XHFHZHX and 


NFBL7 


X^Y^Z^X,Y^X, and Z^Y 


NFBL8 


X^Y^Z^X,Y^X, and X^Z 


NFBL9 


X^Y^Z^X,Z^Y, and X^Z 


NFBLIO 


X^Y^Z^X,X^Y, and Y^Z 



1. Master Equations 

To simplify mathematical analysis, we assume that (a) any given node is constitutively syn- 
thesized, (b) any given node is degraded depending on its current abundance (c) any inhibitory 
interaction is described by mass action kinetics, and (d) the corresponding stochastic chemical 
reaction is a Markovian process. Based on these assumptions, we can write down the chemical 
master equations. To this end, we denote the number of each constituent of the networks X, Y, 
and Z. Thus, we can describe the state of a system as a state vector, i.e., X = {X, Y, Z). We also 
denote the joint probability distribution of a network in a state X at a particular time t as P{X; t) 
and the transition rates between two different states as T{X\X'). Therefore, we can write down 
the master equation governing the time evolution of P{X; t) formally as 

= E nX\X')P{X'; t) - E T{X'\X)P{X; t). (1) 
X' X 

The possible transitions among the states are dependent on the specific reactions allowed in each 
network. In Table [III show how to write down the transition rates of synthesis and degradation 
associated with all the possible chemical reactions that can occur on a single node. 

Based on Table [III can write down all the transition rules for NFBLl as follows. 

T{X + l,Y,Z\X,Y,Z) = kiQ. (2) 
T{X -1,Y,Z\X,Y,Z) = k2XZ/n. 
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TABLE II: The transition rates of synthesis and degradation associated with each chemical 

reaction 



Reactions 


Transition rates for synthesis 


Transition rates for degradation 


X 


kiVL 


k2X 


Y ^ X 


kiY 


k2X 


Y ^X, Z^X 


kifl 


k2XZ/VL 


Y ^X, Z^X 


kifl 


kaiXY + XZ)/Q. 


Y ^ X, Z ^ X 


ki{Y + Z) 


k^X 



T{X,Y + 1,Z\X,Y,Z) = hX. 

T{X,Y -1,Z\X,Y,Z) = kiY. 

T{X,Y,Z + 1\X,Y,Z) = k^Y. 

T{X,Y,Z -l\X,Y,Z) = keZ. 

Since the number of each chemical species can be raised or lowered only by 1 in each reaction 
process, it is convient to introduce the step operator Ej^ and E^, the operation of which is defined 

by 

E^iX, Y, Z) = fix ±1,Y,Z) = Y^ Z) (3) 

m=0 

where f{X, Y, Z) is an analytic function. Then, we can rewrite the master equation by using the 
step operators, 

E [Ex^-mX'\X)PiX)+ [El-mX\X')P{X') (4) 

Xi=X,Y,Z x,=x,Y,z 

where the sum is carried over all possible reactions raising or lowering the number of each species 
in the networks. 

2. Macroscopic rate equation and Fokker-Plank equation 

The master equation involves non-linearity and for such a case, the analytic solution is not 
available in general. Thus, we use the system size expansion method by Van Kampen [35], which 
allows us to expand the master equation systematically in terms of the system size Thereby, we 
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can obtain the approximate solutions to the leading orders of the system size expansion. Following 
the standard procedure discussed in Ref. [35], we define 

Xi = n4>,^ + (5) 

where the index i runs over all the chemical species, i.e., Xi = X, Y, Z. il.(j)xi describes the macro- 
scopic mean value of Xi and the Cl^^'^^i is the gaussian fluctuations of Xi around its mean. Then, 
we can rewrite the joint probability distribution P({Xi};t) in terms of the stochastic variables of 
such that P{{Xi}) = n({,^j}) = n(^j). According to the change of variables, the time derivative 
of the joint probability distribution transforms as 

dPjx,) _ duji) ^ d<i>,^ dnii) 

dt dt ^ dt dii ' ^ ' 

si 

Consequently, the step operators are also rewritten in terms of the new variables. 



Ex^Ef = Y^{±irn-f-^. (7) 

m=0 

Plugging Eq. ^ through Eq. ([7]) to Eq. Q, we can expand the master equation systematically 
in the order of the system size ft. To the leading order of expansion we obtain the 

macroscopic deterministic equations for all the networks. For example, we find that for NFBLl, 
they are given by 

^ = h - k2<l)x(py, (8) 



dt 



k3(t>x - kicjiy, (9) 



= k5(j)y - ke(p^. (10) 

We can simply repeat the same procedure for the rest of the networks to derive the macroscopic 
deterministic equations. In the next leading order of Q, expansion (il*^), we obtain the linear 
Fokker-Plank Equation. 

a^ = -E^^te)n+E2^.ape; (11) 

where C(^j) = Ylj JijCj- Jij denotes the Jacobian matrix of the macroscopic rate equations and 
Bij is the noise covariance matrix. For NFBLl, the stability matrix and the covariance matrix 
evaluated at the fixed points are given by 



J 



k^ —k4 
\ fcs -ke J 



(12) 
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and 



B 







V 



k. 










where 



, and (j)*^ refer to the fixed points of the macroscopic rate equations. 



(13) 



3. Power spectrum 



In order to examine the existence of the stochastic ociUations, it is useful to transform from the 
time domain to the frequency domain, i.e., the Fourier transform. For this purpose, we recast the 
linear Fokker-Plank equation into the mathematically equivalent form of Langevin equations, 

dm ^ 

3 



dt 



(14) 



where < r]i{t)r]j{t') >= Bij5{t — t'). Taking the Fourier transform of Eq. (14), —iuS^i^uj) = 
YjjJijijit) + fiiit) where we introduce ii{uj) = J dte~'''^^ ii{t) and fii{uj) = J dte-^'^^rjiit). We 
then solve for Xi{uj) to get Ci{u}) = ^^j^{uj)rjj{uj) where we define ^ij = —iuSij — Jij. Finally, 
we obtain the power spectrum 



p.H =< Uu)\' >= E^.^-'(^)^.fc(^fc/)^M = ^'^^\u;)B,,{<^j,y{^ 

j,k j 



UJ 



(15) 



where f denotes the conjugate transpose. The last equality comes from Bij being diagonal. Since 



it further simplifies to 



P,(a;) = Y,B,j{<^>-j')i-u;)^-jHu;)B,, = ^ S,,|$-.i(^) 



(16) 



3 J 

For 2 by 2 J and B matrices, the algebraic expression of the power spectrum is simple and according 
to Ref. [31], it is given by 

Ojo;^ + bi 



uj'^ + at<;2 + 13 



where i = 1,2, and 



«i — Bii, bi — B11J22 + B22J12, 



(17) 



(18) 



02 = B: 



22, 



b2 = B11J21 + B22J11, 



and a = Tr[J]'^ - 2Det[J], ^ = Det[Jf 
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In particular, the condition for the existence of power spectrum peak is equivalent to that of 
z > such that dPi{z)/dz = where z = u:'^: aiz'^ + 2biZ + abi — aj/3 = 0. According to the 
Descartes' rule [SB], the condition for the quadratic equation to have a real positive root is given 

by 

abi-aif3<0 ^ {Tr[J]'^ - 2Det[J])b - aDet[jf < (19) 



where Oj > and bi > are trivially satisfied as shown in Eq. (18). 

For our networks, J and B are given by 3 by 3 matrices, the explicit expression of the power 
spectrum for i = x,y,z in terms of matrix elements of J and B is given by 

Fi[uj) = —7. -. 7, (20) 

u}^ + aiu}^ + OL2W + as 

where the coefficients in the numerator are: /3i = Bjj, /32 = Bii\{^r\J\ — Ja)'^ — 2Mjj] + 
+ Bi+2A+2Jli+2, and /Ss = BuM^- + Si+i^j+iM^.^^ + 5^+2,^+2^2.^2 and Mij de- 
notes the minor matrix of Jij. The index of the matrices obeys the following conventions: 
Bij = -Bi+3j+3, and Mij = Mj+sj+a. The coefficients in the denominator are: 

ai = rr[J]2 - 2(JiiJ22 + J22-^33 + -^33-^11 - -^12^21 - J23J32 - J31J13), "2 = {J11J22 + ^22-^33 + 

J33^ii - ^12^21 - ^23^32 " ^31^13)^ " 2Tr[J]Det[J], and as = Det[J]'^. 



B. Linear Stability Analysis 



The set of deterministic rate equations for each network takes the general form of 



(21) 



For such a system, determining the stability of the steady state, 



0, entails the 



examination of the linearised form of (21 ). Namely, in performing a multivariate Taylor expansion 



about the fixed points, {<!)%, (py,(p*}, the above system of non-linear ODEs may be approximated 
by the first order term which is linear with respect to 6(pi = {(pi — (p*). More explicitly, 

/-iri /-lo 



V 



Ju J12 Jl3 
J2I J22 J23 
\ J3I J32 J33 J 



(22) 



V Hz J 



where J,, 



represents the i^^^f^ element of the Jacobian matrix. According to the 



Hartman-Grobman theorem, so long as the eigenvalues of the Jacobian do not lie on the imaginary 
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axis, the linearised form may be used to assess the stabihty of the steady state for the non-Hnear 
system. In light of this fact, we begin then by grouping each feedback loop into a family derived 
from either NFBLl or NFBL2. In doing so, one notices for NFBL2, NFBL6, and NFBLIO, the 
Jacobian takes on the general form 



/ -an ai2 -ais ^ 



J 



(23) 



-O21 — a22 023 

y -032 -033 J 

for ttij = I Jjjl being greater than zero so long as each fixed point is greater than zero. A proof 
of stability for the deterministic system amounts to demonstrating Re{X) < 0, where A denotes 
the eigenvalue of the Jacobian matrix. Finding such eigenvalues requires one to solve for the 
characteristic equation of the 3x3 determinant, det{J — IX) = 0, given here by 

- Tr( J)A2 + (^11 + A22 + ^33)A - det{J) = (24) 

for Tr{J) representing the trace of the Jacobian, while A^j is the cofactor of the z*'* row and 
column. By applying the Routh-Hurwitz conditions for a third order characteristic equation, it is 
guaranteed Re{\) < so long as Tr{J) < 0, det{J) < 0, det{J) - {An + A22 + A2,^)Tr{J) > 0. 

det{J) - {An + A22 + As3)Tr{J) (25) 

= flliai2021 + ai2fl2l022 + 022023^32 + (1230132033 
+ 011(122 + 011022 + OllOsa + 0^2033 

+ 0110^3 + 0220I3 + 2aiia22033 - 013021032. 

Interestingly enough, such an expression will always positive provided 2011022033 — 013O21O32 > 0. 
For each network derived from NFBL2, it turns out 2011022033 — 013O21O32 = k2k4kQ(l)*(j)y(p* for 
any general case. In noting this quantity to always be positive, it follows, that all networks in 
the family of NFBL2 are stable. For the family of the networks stemming from NFBLl, it can 
be shown that all three Ruth-Hurwitz conditions always hold. Consequently, all NFBL networks 
presented have stable steady state solutions. 



C. Numerical methods 



In the following, we explain the details of sampling methods, power spectrum and SNR cal- 
culations. We also explain how we define the robustness and the prominence of the stochastic 
oscillations for our study. 
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1. Sampling methods 

As shown in the stability analysis, our networks are entirely stable for positive real domain 
in the parameter space, i.e., {ki} G In order to demonstrate the existence of the stochastic 
amplified oscillations, we sample 10^ sets of the 5 kinetic rate constants (setting /c2 = 1) uniformly 
and randomly from the logarithmically scaled interval [10^^, 10^] by using the unconstrained Monte 
Carlo sampling technique. We also used three other different sampling methods: MC sampling from 
linear parameter space, LHS (Latin Hypercube sampling) from linear and logarithmic parameter 
space. Both MC and LHS Sampling from linear parameter space didn't provide a good number of 
SNR> 1. The LHS sampling from logarithmic parameter space produced the very similar results 
as the MC sampling from logarithmic parameter space. 



2. Power Spectrum and SNR calculation 

The SNR (signal to noise ratio) is commonly defined by the power spectrum peak height over 
its relative width, i.e., 

SNR = (26) 

where ojo and P{uJo) denote the peak frequency and the peak height respectively. The Aw defines 
the so-called full width at half maximum (FWHM). Our calculation of the power spectrum and 
the SNR is based on its analytic expression obtained from the earlier section. It is, however, 
not straightforward to obtain the analytic expression for SNR because SNR calculation requires 
the solution of 4th order polynomial equation given by dP{z)/dz = where z = cj^ and the 
determination of Ao; around is also non trivial. Therefore, we use a numerical algorithm to 
calculate the SNR in the following way. First, we identify the occurrence of the maxima of P(c<j) for 
a given set of the kinetic rate constants and fixed points. In this step, we obtain the numerical values 



of P{oo) from the analytic expression of Eq. (20) and then examine the condition for the existence 
of maxima. For example, P{oo) values are evaluated at between a; = to 200 in the increment 
of Aw = 0.01. In each increment, we calculate the difference of APi = P{oo) — P{uj — Aoj) 
and AP2 = P{^ + Aw) — P{uj). If APi > and AP2 < 0, we locate the peak frequency at 
uj = (jJo- Second, if the peak frequency Wq exists, we search for the two neighboring frequencies, i.e., 
uJi{< uJo) and W2(> Wq) that satisfy -P(wi) = P(w2) = P{uJo)/2. Then, we can determine Aw from 
the difference of two frequencies, i.e., Aw = |w2 — wi|. 
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3. Probability distribution of SNR 

The probability distributions of SNR is obtained from the histogram of SNR values in the range 
of to 10'' followed by the normalization. To obtain the histogram, we set the bin size to vary in 
the linear scale to smooth out the noisy curves, which occur due to the fluctuations in the number 
of data sets as SNR values get larger. For example, in the range between SNR =1 and 10, the bin 
size is set to 1 and in the range between SNR =10 and 100, the bin size = 10. 

4-. Robustness vs. Prominence 

In this work, the robustness of a networks is defined as the size of the domain of the parameter 
space that generates the stochastic amplified oscillator, i.e., SNR> 1 whereas the prominence is 
the average value of SNR within such domain for that networks. We measure the domain of the 
parameter space from counting the data sets showing SNR > 1 out of total sampling points and 
calculate the average value of SNR within the domain by integrating the probability distribution 
of SNR, i.e., 

{SNR)= [ Pi{SNR)Q{SNR-l)d{SNR) (27) 

i=x,y,z 

where Di denotes the magnitude of the domain that show SNR > 1 for a given species i = x,y, z. 
So Di is equivalent to the counts of the data sets with SNR > 1 for a given species in a model. The 
summation over x, y, z is carried out only among those species that exhibit SNR> 1. For example, 
NFBL3 do not show SNR > 1 for x species. Thus the summation is limited to y and z species. 
The step function 0(x) is inserted to restrict the integral range only for SNR> 1. 

III. RESULTS 

A. All the networks with 3 nodes and mass action kinetics and only negative feedback 
loops are capable of being stochastically amplified oscillators. 

1. Network with two components, a negative feedback, and mass action kinetics 

The networks we consider in this work are built upon the following constraints: i) The networks 
have mass action kinetics ii) No positive feedback loop, iii) 3 nodes. The third constraint, 3 nodes, 
is not a necessary condition for stochastic oscillations as discussed in Ref. [211 132]. But, if we 
add the constraints i) and ii), having 3 nodes appears to be a necessary architectural condition to 
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FIG. 2: Probability distribution of Signal to Noise Ratio (SNR) for (a) a network with 
two-component negative feedback loop and (b)-(d) nine networks with three-component negative 
feedback loop. The SNR distributions for both species X and Y are plotted in (a) whereas those 
for X, Y, and Z species are plotted in (b), (c), and (d), respectively. We note that not all species 

of a network necessarily exhibit SNR>0. 

build a good stochastic oscillator. We can show that NFBL with 2 nodes is not a good stochastic 
oscillator. 

We consider a network with 2 component NFBL constructed with the same design rules used 
for the 3 component NFBL network counterpart: X and Y -\ X. For the detailed discussion 
of the 2 component network and the derivation of its power spectrum, we refer the readers to 
the methods section and Ref. [271. We find that there exist some sets of kinetic rate constants 



that satisfy Eq (19), i.e., the condition for the existence of the peak of the power spectrum. But, 



having a peak does not necessarily mean the existence of highly amplified and coherent stochastic 



14 



• 1— I 

o 

Ph 



I40 

120 
100 

BO 

eo 

40 

20 





1 

! 
i 
i 


Deterministic — 
Stochastic 






















M 



















400 600 

Time 



o 

Oh 
U 

O 





Analytic — 
Numerical o 









Frequency 



FIG. 3: Top panel: comparison between the stochastic time series and the deterministic time 
series for Y species of the NFBL network 1. Bottom panel: comparison between two power 
spectra, one from the linear noise approximation and another from the numerical simulation. 

oscillations. In Fig. [2] (a), we show the probability distribution of SNR for X and Y species for the 
2 component NFBL network. We sample 10^ sets of the kinetic rate constants in the logarithmic 
scale from 10"^ to 10^. The SNR values appear to be bounded below 1. This numerically proves 
that the 2 component NFBL network is not a good model to obtain sufficient amplified stochastic 
oscillations. 

2. Networks with three components, negative feedbacks, and mass action kinetics 

In contrast to the 2 component NFBL network, the 3 component NFBL networks exhibit highly 
amplified and coherent stochastic oscillations. In Fig. [s] (a), we compare a stochastic time series 
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data with the corresponding deterministic trajectory obtained from solving the macroscopic rate 
equations for NFBL network 1. The time series data don't appear to have clearly discernable 
periodic behavior. However, as shown in Fig. [s] (b), we can see that the power spectrum obtained 
from averaging over 1000 realizations of such stochastic time series data shows a prominent peak 
with small width, which is the evidence of the periodicity of stochastic time series data. It also 
verifies that the power spectrum derived from the linear noise approximation is in good agreement 
with that simulated with Gillespies direct algorithm in a carefully chosen systems size. 

As shown in the Fig. [2] (b), (c), and (d), all the networks with 3 component NFBL exhibit 
the stochastic oscillations quantified by the SNR in a broad range from to 10^. All 9 networks 
are deterministically stable in the entire parameter space but noise is capable of inducing highly 
amplified and coherent stochastic oscillations. It appears that the SNR probability distributions 
for all the networks follow a similar trend of a power-law behavior followed by a cut-off. We have 
not investigated if the cut-off is dependent on the sampling size within the biologically relevant 
parameter value range. 

B. Networks with a single negative feedback loop are better stochastic oscillators than 
those with multiple coupled negative feedback loops 

As shown in the previous section A, all networks are capable of being stochastic simplified and 
coherent oscillators if the kinetic rate constants are properly balanced. In this section, we compare 
those 9 networks and identify the network architectural conditions for, namely the design principles 
of, the most prominent and robust stochastic oscillators. 

In Fig.|4]is shown the prominence versus the robustness of all nine networks. For our particular 
networks, the prominence is defined as the average value of the signal to noise ratio (SNR) over 
the stochastic oscillators with SNR> 1 and over three variables {X, Y, and Z) . In this prominence 
measure, we exclude the cases with SNR = and SNR< 1 because the prominence is supposed 
to measure how good the stochastic oscillators are when they satisfy the noise-induced oscillation 
conditions: the existence of a peak of a power spectrum and the signal being greater than stochastic 
fiuctuations, i.e., SNR> 1. The robustness is measured by the average percentage of the cases of 
SNR> 1 over the total sampling size. Each network has three nodes (X, y, and Z) and all three 
variables are not necessarily stochastically oscillatory nor have the same average SNR. For example, 
a variable Y of a network can exhibit a highly amplified stochastic oscillation while a variable X 
of that network cannot have a single case of SNR> 1 among the total samples. In Fig. [4j we 



16 



0.1 



m 

o 
o 



-4— > 
O 



0.01 - 



0.001 



1 

- 
- 

- 
- 


1 1 1 1 ' ' 1 ' 

VVUlaL 


1 1 1 1 1 


1 ' ' ' 1 
Best 


1 1 1 

- 
- 

- 






/ . \ 
1 ) 

\ / 


■ 


X + ) 






^ ^ 




NFBL 1 


+ 
X 








NFBL2 








NFBL 3 










NFBL 4 
NFBL 5 


□ 
■ 








NFBL 6 










NFBL 7 


• 








NFBL 8 


A 


1 


III 1 III 1 


III 1 1 


NFBL 10 
-1 1 1 — 





10 100 1000 10000 100000 19+06 

Prominence (Average SNR) 



FIG. 4: Network design space for identification of stochastically prominent and robust oscillators. 
Robustness is plotted against prominence. Please refer to both main text and methods for the 

detailed definition of robustness and prominence. 

present the prominence and the robustness averaged over the cases of SNR> 1 and over the three 
variables of a network. Please note that two measures are similar to the maximum of the average 
SNR among three variables. Nine networks are classified into three distinctive groups, based upon 
the prominence measure. The NFBL networks 3, 7, and 8 are categorized as the worst performers, 
the NFBL networks 1, 2, and 4 as the best performers, and the NFBL networks 5, 6, and 10 as the 
mediocre performers. Between the best and the worst performers, there is a large difference in the 
prominence measure, by 3 to 4 orders of magnitude. 

We identify the common network architecture among the worst performers: the NFBL networks 
3, 7, and 8 have competing double inhibitions acting on a node X as shown in Fig 1. For all of the 
worst performers, the node X has zero case with SNR> 1 out of the total sample size 10^. When 
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both inhibitions act on the same node in an asynchronous and random manner, they compete 
to influence the node X in different directions, resulting in cancelling out of their driving forces 
on the node. This is similar to the case of a pendulum perturbed by two random forces, whose 
forces are completely cancelled out on average. The similar situation can be found in deterministic 
oscillators. 

All of the worst performing networks belong to the family of the networks derived from NFBL 
network 1 which can be generated by adding one or more directed inhibition or activation edges 
to NFBL network 1 with a single constraint that any additional edge should not create a positive 
feedback loop. Thus, it is topologically allowed in this family of NFBL network 1 to have double 
inhibitions or double activations acting on the same node. However, the family of networks stem- 
ming from NFBL network 2 is in principle not allowed to have any competing double inhibitions 
or activations. Any additional directed edges introduced to the NFBL network 2 are allowed to be 
only an activation. This architectural constraint make this latter family of NFBL networks either 
the best or the mediocre stochastic amplified oscillators. It is notable that the same network ar- 
chitecture as NFBL network 2 was used to demonstrate the synthetically designed gene regulatory 
oscillations in noisy intracellular setting in Escherichia coli |37j . 

C. Multiple timescale difference among the kinetic rate constants is required for stochastic 

amplified and coherent oscillations. 

The kinetics that synthesize and degrade the interacting biochemical species must occur on 
appropriate timescales that permit a network with three-component negative feedback loop to 
generate the stochastically amplified and coherent oscillations. 

In Fig. [5] are shown the two-dimensional phase diagrams for the X species of three best perform- 
ing networks, NFBL networks 1, 2, and 4. The SNR values are averaged over three other kinetic 
rate constants and mapped onto the two-dimensional parameter space. For all three networks, the 
highest average SNR values are located around the left bottom corner of the diagrams where both 
kinetic rate constants are in the order of 10^^ to 10~^. The similar distribution of the average SNR 
values are found, but not presented here, for the X species of three other networks with mediocre 
performance, NFBL networks 5, 6, and 10. As for the Y species, all six networks exhibit the highest 
average SNR values around the right bottom corner where one kinetic rate constant is in the order 
of 10^ and another in the order of 10~'^, thus indicating their ratio in the order of 10^. 

Since we rescale the time with k2 (i.e., setting ^2=1)) the other rate constants are effectively 
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FIG. 5: Phase diagram of Signal to Noise Ratio for X species of three NFBL networks that 
generate the most prominent and robust stochastic oscillations. Those SNR vales are averaged 

over 3 kinetic rates constants and coarse-grained. 

normalized by k2. The characteristic distribution of the average SNR in two dimensional phase 
diagrams in Fig. [5] indicates that it is required for the stochastic amplified and coherent oscillations 
that both kinetic rate constants and /cg for NFBL networks 1 and 4 (or and k^ for NFBL 
networks 2) should be by two to three order of magnitude smaller than the kinetic rate constant 
/c2- Both k^ and k^ are the synthesis rate of Y and Z species and k2 is the degradation rate 
of X species. Thus, this finding suggests that there must be a multiple timescale requirement 
for stochastic oscillations similar to that found in Ref. [H], i.e., proper balancing among rates 
of synthesis and degradation of interacting chemical species for a network with three-component 
negative- feedback loop to generate oscillations. For NFBL network 2, all three rates k2, k^, and fee 
are the degradation rates induced by inhibition. This also hints that in order for NFBL network 
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2 to generate the stochastic oscillations, there must be multiple timescale difference among the 
degradation rates. 

IV. DISCUSSION 

The goal of this paper is to investigate the relationship between network structure and stochastic 
dynamical behaviors of small-size biological networks with coupled negative feedbacks. We identify 
the design principles of the best stochastic amplified and coherent oscillators that are made of 
only negative feedbacks: three components, a single or at most two negative feedbacks, and no 
competing double-inhibitions acting on a node. In other words, the more number of negative 
feedbacks a network has, the more likely it has double-inhibitions, and thus the less likely it 
exhibits stochastic oscillation. This design agrees with the common network architectures found 
among biological oscillators. E.g., the core part of the fungal circadian clock operates with a single 
transcription/translation negative feedback loop where a single gene frq inhibits its own expression 
in a stochastically cyclic manner [6]. There is a counter example, too. Nucleo-cytoplasmic shuttling 
of NF-kB is regulated by multiple coupled negative feedbacks such as IkBo, IkB/3, InBe, and A20. 
This signaling network is not the best stochastic oscillator according to our work, but exhibits 
noisy oscillations in single cells [38) . One possible explanation might be that the NF-kB signaling 
network is coupled with a positive feedback through autocrine signaling of TNFa and thus should 
be regarded as the interconnected positive and negative feedbacks which is one of the most common 
coupled feedbacks in biological oscillators. 

Since any intracellular biochemical reactions are subject to random noise, a plethora of noise- 
induced oscillation phenomena have been found in biological oscillations and extensively investi- 
gated for over 30 years [H [25l - [33l [35} [371 [38] . However, there exist not much theoretical work 
available on the necessary and sufficient conditions for noise- induced oscillations [251 [3T] . In this 
paper, we are intrigued by a similar, but biologically relevant question: what are the necessary net- 
work architectural conditions for noise-induced oscillations? Our numerical work uncovers one of 
those conditions, a single negative feedback. But, we are left with many more questions. What are 
the precise mechanisms that make multiple negative feedback loops detrimental to noise-induced 
oscillation? How can noise induce oscillation when the deterministic dynamical systems have no 
center and thus no limit cycle in the entire parameter space? One extreme case related to the latter 
question is very intriguing. Our unpublished data show that the signal to noise ratio, a metric to 
judge noise-induced oscillation, can be very high even though all of the eigenvalues of the stability 
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matrix are real and negative. In this case, the deterministic system does not operate with any 
natural frequency, but some unknown mechanisms govern and determine a timing (or resonant 
frequency) of noise-induced oscillation. 

In this paper, we consider a small set of networks with three components, negative feedback 
loops, and mass action kinetics. The above three constraints are to simplify and ease our modeling 
exercise and thus reflect the limitations of our models. In our future work, we relax those constraints 
to investigate networks with N components, both positive and negative feedback loops, and/or 
nonlinear kinetics. From our comparative studies between two-component network and three- 
component networks, we can conjecture that the larger number of intermediate components a 
negative feedback network has, the more time-delay it has, and the more likely it exhibits stochastic 
oscillation. This conjecture can be tested with a stochastic version of N-cyclic network with negative 
feedbacks |39l HOJ. Our current modeling framework can be easily modified and extended to 
enumerate networks with interconnected positive and negative feedbacks and identify the most 
robust and prominent stochastic oscillators within this group of networks. 
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